{
"cells": [
{
"cell_type": "markdown",
"id": "c2f4fd6c",
"metadata": {},
"source": [
"# Implementing Shor's algorithm in Perceval"
]
},
{
"cell_type": "markdown",
"id": "5162c9b8",
"metadata": {},
"source": [
"This notebook presents a simulation in Perceval of a 4-qubit 12-modes optical circuit performing Shor's algorithm, based on Alberto Politi, Jonathan C.F. Matthews, and Jeremy L. O'brien. \"Shor’s quantum factoring algorithm on a photonic chip.\" [Science](https://www.science.org/doi/10.1126/science.1173731) 325.5945 (2009): 1221-1221."
]
},
{
"cell_type": "markdown",
"id": "f29e5abf",
"metadata": {},
"source": [
"## Shor's algorithm"
]
},
{
"cell_type": "markdown",
"id": "5dd3a2a5",
"metadata": {},
"source": [
"The purpose of Shor's algorithm is to find nontrivial factors of a given number $N$ in polynomial time, yielding an near-exponential speedup compared to state of the art classical algortihms.\n",
"\n",
"The main routine of Shor's algorithm consists in finding the order $r$ of a number $a \\in \\mathbb{Z}_N$, i.e. the smallest integer $r$ such that $a^r = 1 \\pmod N$.\n",
"\n",
"If the order of a randomly chosen $a$ which is coprime with $N$ is even, then $(a^{r/2} - 1)(a^{r/2} + 1) = k N$. If none of these factors are multiples of $N$, then $\\gcd(N, a^{r/2} - 1)$ and $\\gcd(N, a^{r/2} + 1)$ are nontrivial factors of $N$."
]
},
{
"cell_type": "markdown",
"id": "ff756302",
"metadata": {},
"source": [
"## Preliminaries"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "55fe7a18",
"metadata": {},
"outputs": [],
"source": [
"import sympy as sp\n",
"import perceval as pcvl"
]
},
{
"cell_type": "markdown",
"id": "826377de",
"metadata": {},
"source": [
"### Path encoding functions\n",
"\n",
"The following functions allow for conversion between the qubit and Fock state representations."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "cd2bbda8",
"metadata": {},
"outputs": [],
"source": [
"def toFockState(qubitState):\n",
" # path encoding\n",
" pe = {0:[1,0], 1:[0,1]}\n",
" return [0] + pe[qubitState[0]] + pe[qubitState[2]] + [0, 0] + pe[qubitState[1]] + pe[qubitState[3]] + [0]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "55ee4ce8",
"metadata": {},
"outputs": [],
"source": [
"def toQubitState(fockState):\n",
" # qubit modes\n",
" x1 = [1, 2]\n",
" f1 = [3, 4]\n",
" x2 = [7, 8]\n",
" f2 = [9, 10]\n",
" # auxiliary modes\n",
" am1 = [0, 5]\n",
" am2 = [6, 11]\n",
" \n",
" # auxiliary modes\n",
" for i in am1 + am2:\n",
" if fockState[i]!=0:\n",
" return None\n",
" L=[]\n",
" # qubit modes\n",
" for q in [x1, x2, f1, f2]:\n",
" if fockState[q[0]]+fockState[q[1]] != 1:\n",
" return None\n",
" else:\n",
" L.append(fockState[q[1]])\n",
" return L"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "39b1c269",
"metadata": {},
"outputs": [],
"source": [
"def strState(state):\n",
" return str(pcvl.BasicState(state))"
]
},
{
"cell_type": "markdown",
"id": "016e2534",
"metadata": {},
"source": [
"## The circuit"
]
},
{
"cell_type": "markdown",
"id": "01c60397",
"metadata": {},
"source": [
"### Quantum circuit\n",
"\n",
"The quantum circuit has been optimized after choosing parameters $N = 15$ and $a = 2$, and aims to calculate $r=4$.\n",
"It features 5 qubits labelled $x_0, x_1, x_2$ and $f_1, f_2$. Qubits $x_i$ and $f_1$ are initially in state $|0\\rangle$, and $f_2$ in state $|1\\rangle$.\n",
"In the non-optimised Shor algorithm, qubits $x_i$ encode a binary number representing a pre-image of the Modular Exponentiation Function (MEF) $x \\mapsto a^x \\pmod N$, while qubits $f_i$ hold the image obtained after applying the MEF to qubits $x_i$. Applying the MEF when qubits $x_i$ hold a superposition of different pre-images (obtained with H gates on qubits $x_i$) allows to efficiently compute the order $r$ of parameter $a$ modulo $N$.\n",
"\n",
"The circuit consists of $\\mathsf{H}$ gates being first applied to each $x_i$ qubit, followed by $\\mathsf{CNOT}$ gates applied on $x_1, f_1$ and $x_2, f_2$ pairs, where $x_i$ are control qubits; finally the inverse $\\mathsf{QFT}$ algorithm is applied on qubits $x_i$.\n",
"\n",
"$\\mathsf{CNOT}$ gates on $x_i, f_i$ pairs ($x_i$ being the control) are implemented using $\\mathsf{H}$ and $\\mathsf{CZ}$ gates: the $\\mathsf{CZ}$ gate is sandwiched between two applications of $\\mathsf{H}$ on $f_i$.\n",
"\n",
"The input state of the circuit after optimisation is $|0\\rangle_{x_0}|0\\rangle_{x_1}|0\\rangle_{x_2}|0\\rangle_{f_1}|1\\rangle_{f_2}$.\n",
"\n",
"The expected output state is then $\\frac{1}{2} |0\\rangle_{x_0} \\otimes \\left ( |0\\rangle_{x_1}|0\\rangle_{f_1} + |1\\rangle_{x_1}|1\\rangle_{f_1} \\right ) \\otimes \\left ( |0\\rangle_{x_2}|1\\rangle_{f_2} + |1\\rangle_{x_2}|0\\rangle_{f_2} \\right )$."
]
},
{
"cell_type": "markdown",
"id": "b7b25cf5",
"metadata": {},
"source": [
"### Photonic circuit\n",
"\n",
"The optical circuit from the result by Politi et al features twelve modes (ordered from 0 to 11 from top to bottom).\n",
"\n",
"During the execution, qubit $x_0$ remains unentangled from the other qubits. It can therefore be removed from the optical implementation.\n",
"\n",
"The qubits $x_1, x_2, f_1, f_2$ are path encoded as modes $(1, 2)$, $(3, 4)$, $(7, 8)$, $(9, 10)$ respectively. The four remaining modes are used as auxiliary modes to implement the $\\mathsf{CZ}$ gates.\n",
"\n",
"With path encoding each $\\mathsf{H}$ gate in the quantum circuit is implemented with a beam splitter with reflectivity $R=1/2$ between the two pathes corresponding to the qubit. In our implementation in Perceval, phase shifters are added to properly tune the phase between each path.\n",
"\n",
"$\\mathsf{CZ}$ gates are implemented with three beam splitters with reflectivity $R=2/3$ acting on six modes: one inner BS creates interference between the two qubits, and two outer BS balance detection probability using auxiliary modes.\n",
"This optical implementation succesfully yields the output state produced by a $\\mathsf{CZ}$ gate with probability 1/9; otherwise it creates a dummy state, which can be removed by post-selection.\n",
"\n",
"In the case $r=4$ the QFT can be performed classically and doesn't need to be implemented in the photonic circuit."
]
},
{
"cell_type": "markdown",
"id": "3acb2273",
"metadata": {},
"source": [
"## In perceval"
]
},
{
"cell_type": "markdown",
"id": "7aa60c25",
"metadata": {},
"source": [
"### Implementing the circuit"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "2fe3a053",
"metadata": {},
"outputs": [],
"source": [
"circ = pcvl.Circuit(12)\n",
"\n",
"# qubit modes\n",
"# for qubit states 0, 1\n",
"x1 = [1, 2]\n",
"f1 = [3, 4]\n",
"x2 = [7, 8]\n",
"f2 = [9, 10]\n",
"# auxiliary modes\n",
"am1 = [0, 5]\n",
"am2 = [6, 11]\n",
"\n",
"\n",
"# H gates\n",
"for q in [x1, f1, x2, f2]:\n",
" circ.add(q, pcvl.BS.H())\n",
"\n",
"# CZ gates\n",
"for x, f, am in [(x1, f1, am1), (x2, f2, am2)]:\n",
" circ.add((am[0], x[0]), pcvl.BS(pcvl.BS.r_to_theta(1/3))) # R = 1/3\n",
" circ.add((x[1], f[0]), pcvl.BS(pcvl.BS.r_to_theta(1/3)))\n",
" circ.add((f[1], am[1]), pcvl.BS(pcvl.BS.r_to_theta(1/3)))\n",
"\n",
"# H gates\n",
"for q in [f1, f2]:\n",
" circ.add(q, pcvl.BS.H())"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "034d5e45",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optical circuit for Shor's algorithm\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
""
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(\"Optical circuit for Shor's algorithm\")\n",
"pcvl.pdisplay(circ)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "00d193de",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The associated matrix\n"
]
},
{
"data": {
"text/html": [
"$\\left[\\begin{array}{cccccccccccc}\\frac{\\sqrt{3}}{3} & \\frac{\\sqrt{3} i}{3} & \\frac{\\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\\frac{\\sqrt{6} i}{3} & \\frac{\\sqrt{6}}{6} & \\frac{\\sqrt{6}}{6} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & \\frac{\\sqrt{6}}{6} & - \\frac{\\sqrt{6}}{6} & \\frac{\\sqrt{3} i}{3} & \\frac{\\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & \\frac{\\sqrt{6} i}{6} & - \\frac{\\sqrt{6} i}{6} & \\frac{\\sqrt{3}}{3} & 0 & \\frac{\\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & \\frac{\\sqrt{6} i}{6} & - \\frac{\\sqrt{6} i}{6} & 0 & \\frac{\\sqrt{3}}{3} & - \\frac{\\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & \\frac{\\sqrt{3} i}{3} & - \\frac{\\sqrt{3} i}{3} & \\frac{\\sqrt{3}}{3} & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & \\frac{\\sqrt{3}}{3} & \\frac{\\sqrt{3} i}{3} & \\frac{\\sqrt{3} i}{3} & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & \\frac{\\sqrt{6} i}{3} & \\frac{\\sqrt{6}}{6} & \\frac{\\sqrt{6}}{6} & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \\frac{\\sqrt{6}}{6} & - \\frac{\\sqrt{6}}{6} & \\frac{\\sqrt{3} i}{3} & \\frac{\\sqrt{3} i}{3} & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \\frac{\\sqrt{6} i}{6} & - \\frac{\\sqrt{6} i}{6} & \\frac{\\sqrt{3}}{3} & 0 & \\frac{\\sqrt{3} i}{3}\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \\frac{\\sqrt{6} i}{6} & - \\frac{\\sqrt{6} i}{6} & 0 & \\frac{\\sqrt{3}}{3} & - \\frac{\\sqrt{3} i}{3}\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\frac{\\sqrt{3} i}{3} & - \\frac{\\sqrt{3} i}{3} & \\frac{\\sqrt{3}}{3}\\end{array}\\right]$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"print(\"The associated matrix\")\n",
"pcvl.pdisplay(circ.compute_unitary(use_symbolic=False))"
]
},
{
"cell_type": "markdown",
"id": "7bdc227b",
"metadata": {},
"source": [
"### Input state"
]
},
{
"cell_type": "markdown",
"id": "4e1ced97",
"metadata": {},
"source": [
"In the preliminaries we provided functions to map qubit states to the corresponding Fock states and vice-versa.\n",
"\n",
"A *computational basis qubit state* on qubits $x_1, f_1, x_2, f_2$ is represented with a list of 4 boolean values.\n",
"\n",
"A *Fock state* is represented with a list of twelve integer values.\n",
"As described above, for Fock states, the modes are enumerated as follows:\n",
"* mode pairs $(1,2), (3,4), (7,8), (9,10)$ respectively correspond to states $0,1$ for qubits $x_1, f_1, x_2, f_2$\n",
"* modes $0,5,6,11$ are auxiliary modes used for CZ gates"
]
},
{
"cell_type": "markdown",
"id": "51d9eca7",
"metadata": {},
"source": [
"The input state of the circuit is $|0\\rangle_{x_1}|0\\rangle_{x_2}|0\\rangle_{f_1}|1\\rangle_{f_2}$.\n",
"With path encoding this corresponds to sending 4 photons in total in the optical circuit, in the input modes corresponding to the inital state of each qubit."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "79f7ca45",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Input qubit state: |0,0,0,1>\n",
"Corresponding input Fock state: |0,1,0,1,0,0,0,1,0,0,1,0>\n"
]
}
],
"source": [
"qubit_istate = [0,0,0,1]\n",
"istate = toFockState(qubit_istate)\n",
"\n",
"print(\"Input qubit state:\", strState(qubit_istate))\n",
"print(\"Corresponding input Fock state:\", strState(istate))"
]
},
{
"cell_type": "markdown",
"id": "d3f88c02",
"metadata": {},
"source": [
"## Simulation"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "ca10ff5e",
"metadata": {},
"outputs": [],
"source": [
"backend = pcvl.BackendFactory().get_backend(\"Naive\")\n",
"backend.set_circuit(circ)\n",
"backend.set_input_state(pcvl.BasicState(istate))"
]
},
{
"cell_type": "markdown",
"id": "e81a4f44",
"metadata": {},
"source": [
"### Computing the output state"
]
},
{
"cell_type": "markdown",
"id": "c4e62d76",
"metadata": {},
"source": [
"Using perceval we compute the output state obtained with the optical circuit."
]
},
{
"cell_type": "markdown",
"id": "2f3e6f17",
"metadata": {},
"source": [
"#### Amplitudes\n",
"\n",
"We first decompose the output state in the computational basis and plot the corresponding amplitudes."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "0cd16853",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Output state amplitudes: (post-selected on qubit states, not renormalized)\n",
"|x1,x2,f1,f2>\n",
"|0,0,0,0> 0j\n",
"|0,0,0,1> (0.055555555555555566+0j)\n",
"|0,0,1,0> 0j\n",
"|0,0,1,1> 0j\n",
"|0,1,0,0> (-0.05555555555555557-1.734723475976807e-18j)\n",
"|0,1,0,1> (-1.3877787807814457e-17-1.734723475976807e-18j)\n",
"|0,1,1,0> 0j\n",
"|0,1,1,1> 0j\n",
"|1,0,0,0> 0j\n",
"|1,0,0,1> (-1.3877787807814457e-17+0j)\n",
"|1,0,1,0> 0j\n",
"|1,0,1,1> (-0.05555555555555557-1.734723475976807e-18j)\n",
"|1,1,0,0> (1.3877787807814457e-17+0j)\n",
"|1,1,0,1> 0j\n",
"|1,1,1,0> (0.05555555555555559+0j)\n",
"|1,1,1,1> (1.3877787807814457e-17+5.204170427930421e-18j)\n"
]
}
],
"source": [
"output_qubit_states = [\n",
" [x1,x2,f1,f2]\n",
" for x1 in [0,1] for x2 in [0,1] for f1 in [0,1] for f2 in [0,1]\n",
"]\n",
"\n",
"print(\"Output state amplitudes: (post-selected on qubit states, not renormalized)\")\n",
"print(\"|x1,x2,f1,f2>\")\n",
"for oqstate in output_qubit_states:\n",
" ostate = toFockState(oqstate)\n",
" a = backend.prob_amplitude(pcvl.BasicState(ostate))\n",
" print(strState(oqstate), a)"
]
},
{
"cell_type": "markdown",
"id": "43e94f97",
"metadata": {},
"source": [
"The amplitudes obtained with perceval correspond to the expected output state\n",
"$$\\frac{1}{2} \\left ( |0\\rangle_{x_1}|0\\rangle_{f_1} + |1\\rangle_{x_1}|1\\rangle_{f_1} \\right ) \\otimes \\left ( |0\\rangle_{x_2}|1\\rangle_{f_2} + |1\\rangle_{x_2}|0\\rangle_{f_2} \\right )$$\n",
"up to numerical precision, without renormalization."
]
},
{
"cell_type": "markdown",
"id": "300664ac",
"metadata": {},
"source": [
"#### Distribution\n",
"\n",
"We now compute the probabilities to obtain each outcome corresponding to measuring the expected output state in the computational basis."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "2277ec12",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Output state distribution: (post-selected on expected qubit states, not renormalized)\n",
"|x1,x2,f1,f2>\n"
]
},
{
"data": {
"text/html": [
"
\n",
"\n",
"
|0,0,0,1>
|0,1,0,0>
|1,0,1,1>
|1,1,1,0>
\n",
"\n",
"\n",
"
|0,0,0,1>
0.003086
0.003086
0.003086
0.003086
\n",
"\n",
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"input_states = {\n",
" pcvl.BasicState(pcvl.BasicState(istate)): strState(qubit_istate)\n",
"}\n",
"\n",
"expected_output_states = {\n",
" pcvl.BasicState(toFockState([x1,x2,x1,1-x2])): strState([x1,x2,x1,1-x2])\n",
" for x1 in [0,1] for x2 in [0,1]\n",
"}\n",
"\n",
"p = pcvl.Processor(\"Naive\", circ)\n",
"\n",
"ca = pcvl.algorithm.Analyzer(p, input_states, expected_output_states)\n",
"ca.compute()\n",
"\n",
"print(\"Output state distribution: (post-selected on expected qubit states, not renormalized)\")\n",
"print(\"|x1,x2,f1,f2>\")\n",
"pcvl.pdisplay(ca)"
]
},
{
"cell_type": "markdown",
"id": "5e464e2c",
"metadata": {},
"source": [
"The distribution computed with Perceval is uniform over each outcome, which corresponds to the expected distribution obtained in the paper when $x_0 = 0$."
]
},
{
"cell_type": "markdown",
"id": "a79723d8",
"metadata": {},
"source": [
"### Interpretation of the outcomes\n",
"\n",
"For each outcome we consider the values of qubits $x_2, x_1, x_0$ (where $x_0$ is 0) which represent a binary number between 0 and 7, here corresponding to 0, 4, 2 and 6 in the order of the table above.\n",
"After sampling the circuit, obtaining outcomes 2 or 6 allows to successfully compute the order $r=4$.\n",
"Obtaining outcome 0 is an expected failure of the quantum circuit, inherent to Shor's algorithm.\n",
"Outcome 4 is an expected failure as well, as it only allows to compute the trivial factors 1 and 15.\n",
"\n",
"Since the distribution is uniform the circuit successfully yields a successful outcome with probability 1/2. This probability can be amplified exponentially close to $1$ by sampling the circuit multiple times."
]
},
{
"cell_type": "markdown",
"id": "7928b50b",
"metadata": {},
"source": [
"## Reference\n",
"\n",
"> Alberto Politi, Jonathan C.F. Matthews, and Jeremy L. O'brien. \"Shor’s quantum factoring algorithm on a photonic chip.\" [Science](https://www.science.org/doi/10.1126/science.1173731) 325.5945 (2009): 1221-1221."
]
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}